options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

explanatory variable obs.x have noise

xN(x0.sx),yN(b0+b1*x0,s)

ex8-0.stan

normal regression without explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

normal regression with explanatory variable’s noise

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -428.006 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       22      -29.9459    8.1469e-05    0.00165913     0.02854           1       25    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -29.95
##    b0         5.36
##    b1         1.87
##    s          2.71
##    m0[1]     38.35
##    m0[2]     39.07
##    m0[3]     22.51
##    m0[4]     28.55
##    m0[5]     38.94
##    m0[6]     13.39
##    m0[7]     29.06
##    m0[8]      7.33
##    m0[9]     10.49
##    m0[10]    22.58
##    m0[11]    10.50
##    m0[12]    11.12
##    m0[13]    32.81
##    m0[14]    10.93
##    m0[15]    22.21
##    m0[16]     8.77
##    m0[17]    19.08
##    m0[18]    10.43
##    m0[19]    42.56
##    m0[20]    38.11
##    y0[1]     38.88
##    y0[2]     43.32
##    y0[3]     25.03
##    y0[4]     33.24
##    y0[5]     42.69
##    y0[6]     14.89
##    y0[7]     30.94
##    y0[8]      6.32
##    y0[9]     12.85
##    y0[10]    20.92
##    y0[11]    12.45
##    y0[12]    14.17
##    y0[13]    32.38
##    y0[14]    12.64
##    y0[15]    22.58
##    y0[16]     3.00
##    y0[17]    20.43
##    y0[18]     8.41
##    y0[19]    38.89
##    y0[20]    36.65
##    m1[1]      7.33
##    m1[2]     11.25
##    m1[3]     15.16
##    m1[4]     19.08
##    m1[5]     22.99
##    m1[6]     26.91
##    m1[7]     30.82
##    m1[8]     34.73
##    m1[9]     38.65
##    m1[10]    42.56
##    y1[1]      8.33
##    y1[2]     12.01
##    y1[3]     14.24
##    y1[4]     15.40
##    y1[5]     26.49
##    y1[6]     27.86
##    y1[7]     32.41
##    y1[8]     38.91
##    y1[9]     36.88
##    y1[10]    49.14
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -30.47 -30.15 1.23 1.03 -32.82 -29.11 1.01      536      808
##    b0       5.33   5.34 1.13 1.12   3.47   7.12 1.01      294      386
##    b1       1.87   1.87 0.10 0.10   1.71   2.05 1.00      467      975
##    s        3.08   3.02 0.55 0.52   2.32   4.09 1.00     1159      984
##    m0[1]   38.34  38.33 1.14 1.10  36.43  40.23 1.00     1852     1199
##    m0[2]   39.06  39.05 1.17 1.13  37.11  41.03 1.00     1777     1246
##    m0[3]   22.48  22.50 0.69 0.67  21.35  23.59 1.00      750      970
##    m0[4]   28.53  28.53 0.78 0.72  27.24  29.82 1.00     1977     1081
##    m0[5]   38.93  38.92 1.17 1.12  36.99  40.89 1.00     1790     1230
##    m0[6]   13.36  13.36 0.83 0.80  12.03  14.70 1.01      338      491
##    m0[7]   29.04  29.04 0.79 0.73  27.73  30.35 1.00     2101     1077
##    m0[8]    7.30   7.31 1.05 1.06   5.59   8.98 1.01      296      389
##    m0[9]   10.46  10.46 0.92 0.94   9.00  11.94 1.01      313      388
##    m0[10]  22.56  22.58 0.69 0.67  21.43  23.67 1.00      759      970
##    m0[11]  10.47  10.47 0.92 0.93   9.01  11.95 1.01      313      386
##    m0[12]  11.09  11.09 0.90 0.91   9.66  12.54 1.01      319      413
##    m0[13]  32.79  32.80 0.91 0.87  31.27  34.31 1.00     2433     1076
##    m0[14]  10.90  10.91 0.91 0.91   9.46  12.36 1.01      318      419
##    m0[15]  22.18  22.19 0.69 0.67  21.03  23.29 1.00      718      970
##    m0[16]   8.74   8.75 0.99 1.00   7.14  10.34 1.01      303      405
##    m0[17]  19.06  19.06 0.70 0.67  17.93  20.19 1.00      489      777
##    m0[18]  10.40  10.40 0.93 0.94   8.93  11.88 1.01      313      388
##    m0[19]  42.55  42.54 1.33 1.30  40.32  44.82 1.00     1475     1295
##    m0[20]  38.09  38.09 1.13 1.09  36.20  39.97 1.00     1877     1155
##    y0[1]   38.30  38.33 3.30 3.24  32.84  43.56 1.00     1968     1895
##    y0[2]   39.01  38.91 3.30 3.03  33.62  44.47 1.00     1591     1758
##    y0[3]   22.56  22.58 3.26 3.09  17.22  27.94 1.00     2023     2013
##    y0[4]   28.57  28.48 3.26 3.10  23.31  33.81 1.00     1958     2013
##    y0[5]   38.91  38.81 3.34 3.16  33.40  44.36 1.00     2056     1904
##    y0[6]   13.44  13.49 3.27 3.15   8.00  18.68 1.00     1251     1799
##    y0[7]   29.13  29.08 3.20 3.05  23.89  34.35 1.00     1990     1901
##    y0[8]    7.38   7.40 3.25 3.03   2.05  12.68 1.00     1473     1680
##    y0[9]   10.46  10.50 3.20 3.05   5.30  15.64 1.00     1224     1623
##    y0[10]  22.55  22.53 3.28 3.11  17.20  28.16 1.00     1687     2055
##    y0[11]  10.50  10.51 3.34 3.38   4.72  15.90 1.00     1549     1808
##    y0[12]  11.05  10.89 3.22 3.13   6.15  16.33 1.00     1455     1891
##    y0[13]  32.71  32.76 3.20 2.98  27.57  38.02 1.00     2125     1959
##    y0[14]  10.86  10.85 3.22 3.12   5.66  16.17 1.00     1519     1924
##    y0[15]  22.21  22.16 3.27 3.04  16.92  27.79 1.00     1992     1835
##    y0[16]   8.79   8.83 3.30 3.02   3.36  14.25 1.00     1263     1770
##    y0[17]  19.12  19.12 3.26 3.12  13.65  24.39 1.00     1942     1974
##    y0[18]  10.27  10.39 3.16 2.91   5.07  15.35 1.00     1699     1815
##    y0[19]  42.53  42.54 3.36 3.06  37.09  48.02 1.00     1817     1838
##    y0[20]  38.04  38.11 3.28 3.20  32.74  43.38 1.00     2046     2009
##    m1[1]    7.30   7.31 1.05 1.06   5.59   8.98 1.01      296      389
##    m1[2]   11.22  11.22 0.90 0.90   9.80  12.65 1.01      320      424
##    m1[3]   15.13  15.13 0.78 0.76  13.88  16.40 1.01      365      563
##    m1[4]   19.05  19.05 0.70 0.67  17.92  20.19 1.00      488      777
##    m1[5]   22.97  22.99 0.69 0.67  21.83  24.10 1.00      808     1002
##    m1[6]   26.88  26.90 0.74 0.69  25.65  28.10 1.00     1551     1057
##    m1[7]   30.80  30.81 0.85 0.80  29.41  32.22 1.00     2414      996
##    m1[8]   34.72  34.73 0.99 0.94  33.07  36.38 1.00     2265     1127
##    m1[9]   38.64  38.63 1.15 1.12  36.71  40.56 1.00     1821     1198
##    m1[10]  42.55  42.54 1.33 1.30  40.32  44.82 1.00     1475     1295
##    y1[1]    7.33   7.31 3.22 3.10   2.02  12.66 1.00     1518     1746
##    y1[2]   11.18  11.20 3.14 3.09   6.10  16.24 1.00     1415     1732
##    y1[3]   15.14  15.01 3.12 2.96  10.09  20.35 1.00     1687     1655
##    y1[4]   19.12  19.12 3.34 3.12  13.70  24.73 1.00     1965     1892
##    y1[5]   23.00  23.11 3.24 3.04  17.81  28.23 1.00     1915     1845
##    y1[6]   26.88  26.81 3.19 3.03  21.69  32.12 1.00     1950     1701
##    y1[7]   30.81  30.81 3.21 3.00  25.52  36.14 1.00     2220     2033
##    y1[8]   34.72  34.75 3.22 3.14  29.61  39.86 1.00     2009     1930
##    y1[9]   38.54  38.48 3.38 3.19  32.87  43.77 1.00     2120     1991
##    y1[10]  42.35  42.47 3.31 3.18  36.77  47.51 1.00     2067     1925

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 finished in 0.5 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 0.6 seconds.
## Chain 4 finished in 0.5 seconds.
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.8 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.6 seconds.
## Total execution time: 0.8 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd   mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -26.37 -28.99 12.60 10.26 -42.98  2.70 1.06       66       34
##    b0       5.29   5.32  1.36  1.26   2.89  7.48 1.03      193       37
##    b1       1.87   1.87  0.12  0.11   1.67  2.09 1.03      180       33
##    s        2.20   2.29  0.97  1.06   0.59  3.72 1.04       79      147
##    sx       1.01   1.03  0.60  0.71   0.11  1.99 1.07       38       35
##    x0[1]   18.15  17.98  0.85  0.70  16.99 19.78 1.03      243       97
##    x0[2]   17.40  17.52  0.86  0.88  15.91 18.56 1.02      264     1265
##    x0[3]    8.46   8.56  0.85  0.90   7.01  9.61 1.02      211      755
##    x0[4]   12.28  12.35  0.64  0.52  11.17 13.25 1.03     1387     2321
##    x0[5]   17.45  17.58  0.82  0.75  16.00 18.64 1.01      425     1486
##    x0[6]    3.25   3.42  1.09  1.24   1.25  4.60 1.04       89      142
##    x0[7]   13.99  13.77  1.26  1.49  12.43 16.22 1.04       86      251
##    x0[8]    0.90   0.98  0.72  0.56  -0.40  2.02 1.02      731     1163
##    x0[9]    3.18   3.05  0.80  0.66   2.05  4.60 1.01      517      428
##    x0[10]   9.44   9.38  0.70  0.57   8.35 10.60 1.01      820     1865
##    x0[11]   1.56   1.75  1.18  1.36  -0.54  3.02 1.04       83      116
##    x0[12]   3.98   3.84  0.95  1.10   2.73  5.66 1.02      216      526
##    x0[13]  14.84  14.76  0.70  0.56  13.72 16.02 1.02      818     1283
##    x0[14]   3.43   3.29  0.78  0.68   2.33  4.78 1.01      421     1578
##    x0[15]   9.89   9.76  0.96  1.07   8.63 11.53 1.03      141      558
##    x0[16]   1.63   1.72  0.71  0.58   0.38  2.70 1.01      596     1370
##    x0[17]   7.98   7.87  0.81  0.85   6.85  9.35 1.02      278      758
##    x0[18]   2.34   2.48  0.77  0.64   0.96  3.43 1.02      398      591
##    x0[19]  20.01  19.95  0.79  0.58  18.78 21.48 1.05      485      177
##    x0[20]  17.31  17.39  0.73  0.61  16.07 18.49 1.02      623      526
##    m0[1]   39.18  39.34  1.52  1.48  36.57 41.54 1.00      667     2460
##    m0[2]   37.80  37.68  1.77  1.91  35.19 40.71 1.03       87      638
##    m0[3]   21.13  21.17  1.60  1.78  18.64 23.61 1.02      278     1790
##    m0[4]   28.25  28.30  1.22  1.05  26.20 30.21 1.01     2298     2079
##    m0[5]   37.91  37.84  1.69  1.82  35.42 40.71 1.02      137      718
##    m0[6]   11.40  11.40  1.96  2.38   8.42 14.36 1.02      223     1917
##    m0[7]   31.41  31.25  2.25  2.91  28.17 34.82 1.03      131     1101
##    m0[8]    7.00   7.01  1.44  1.29   4.49  9.33 1.02      472       41
##    m0[9]   11.25  11.33  1.55  1.56   8.58 13.56 1.02      192       44
##    m0[10]  22.93  22.91  1.31  1.17  20.90 25.04 1.01      764     1911
##    m0[11]   8.26   8.25  2.11  2.55   5.10 11.58 1.02      216     1389
##    m0[12]  12.74  12.79  1.92  2.24   9.45 15.56 1.04       87       36
##    m0[13]  33.01  33.05  1.28  1.18  30.88 35.08 1.01     2718     2539
##    m0[14]  11.72  11.79  1.57  1.50   8.92 14.16 1.02      155       37
##    m0[15]  23.77  23.72  1.78  2.14  21.16 26.46 1.03      115     1771
##    m0[16]   8.37   8.36  1.39  1.24   6.12 10.66 1.02      502      115
##    m0[17]  20.20  20.21  1.57  1.73  17.62 22.63 1.02      150       91
##    m0[18]   9.69   9.64  1.42  1.35   7.53 12.01 1.01      485     1817
##    m0[19]  42.66  42.68  1.50  1.37  40.16 45.12 1.01      840      729
##    m0[20]  37.63  37.57  1.43  1.37  35.33 40.00 1.01      563     2365
##    y0[1]   39.18  39.45  2.77  2.32  34.24 43.30 1.00     2366     3448
##    y0[2]   37.83  37.40  3.04  2.69  33.54 43.34 1.01      547     1345
##    y0[3]   21.13  20.76  2.83  2.54  17.10 26.12 1.01      875     2472
##    y0[4]   28.26  28.20  2.72  2.24  23.97 32.84 1.01     4028     3276
##    y0[5]   37.91  37.58  2.97  2.60  33.48 43.28 1.01      675     2143
##    y0[6]   11.41  10.96  3.10  2.97   7.17 17.10 1.01      505     2636
##    y0[7]   31.36  31.83  3.34  3.28  25.37 35.84 1.02      334     1723
##    y0[8]    7.00   6.94  2.82  2.33   2.32 11.64 1.01     1439      698
##    y0[9]   11.25  11.51  2.86  2.50   6.22 15.52 1.01     1077     1096
##    y0[10]  22.90  23.07  2.70  2.25  18.35 27.17 1.01     4267     3342
##    y0[11]   8.27   7.77  3.25  3.11   3.75 14.11 1.01      512     2299
##    y0[12]  12.82  13.28  3.04  2.71   7.10 17.08 1.01      358     1744
##    y0[13]  32.98  32.99  2.71  2.31  28.51 37.38 1.01     4057     3615
##    y0[14]  11.76  12.06  2.85  2.43   6.71 16.03 1.01      896     1180
##    y0[15]  23.73  24.17  3.01  2.78  18.24 28.04 1.01      502     2187
##    y0[16]   8.37   8.27  2.88  2.30   3.65 13.31 1.01     1775     1309
##    y0[17]  20.21  20.53  2.85  2.46  15.16 24.39 1.01      750     2128
##    y0[18]   9.68   9.43  2.79  2.37   5.39 14.49 1.01     2317     2746
##    y0[19]  42.66  42.65  2.83  2.39  38.01 47.45 1.01     3075     2787
##    y0[20]  37.55  37.38  2.76  2.31  33.11 42.31 1.00     2921     2568
##    m1[1]    7.26   7.29  1.26  1.17   5.03  9.26 1.03      198       37
##    m1[2]   11.18  11.19  1.06  0.99   9.32 12.84 1.02      215       36
##    m1[3]   15.10  15.11  0.89  0.86  13.53 16.49 1.02      256       35
##    m1[4]   19.02  19.02  0.77  0.76  17.70 20.26 1.01      398       81
##    m1[5]   22.93  22.92  0.72  0.70  21.81 24.14 1.01      723     1035
##    m1[6]   26.85  26.83  0.76  0.74  25.65 28.10 1.01      872      976
##    m1[7]   30.77  30.78  0.87  0.85  29.37 32.18 1.01      672      773
##    m1[8]   34.69  34.70  1.03  1.02  32.99 36.31 1.02      457      717
##    m1[9]   38.61  38.61  1.23  1.18  36.55 40.61 1.02      343      282
##    m1[10]  42.53  42.52  1.44  1.38  40.10 44.85 1.02      274      218
##    x1[1]    1.06   1.06  1.17  0.77  -0.90  2.99 1.03     3896     2269
##    x1[2]    3.13   3.13  1.20  0.76   1.07  5.12 1.03     3995     2078
##    x1[3]    5.25   5.25  1.15  0.79   3.25  7.15 1.03     3794     2527
##    x1[4]    7.33   7.34  1.16  0.74   5.41  9.28 1.03     4198     2082
##    x1[5]    9.42   9.43  1.19  0.81   7.37 11.39 1.03     3848     2244
##    x1[6]   11.52  11.54  1.14  0.75   9.62 13.38 1.03     3596     1975
##    x1[7]   13.60  13.62  1.18  0.77  11.61 15.55 1.03     4012     1802
##    x1[8]   15.73  15.73  1.17  0.76  13.74 17.71 1.03     3716     1926
##    x1[9]   17.83  17.83  1.18  0.80  15.83 19.81 1.03     4280     2634
##    x1[10]  19.92  19.92  1.17  0.78  17.94 21.93 1.03     3803     2127
##    y1[1]    7.22   7.24  2.69  2.33   2.76 11.57 1.00     1187      655
##    y1[2]   11.17  11.16  2.60  2.28   6.94 15.34 1.01     1498     1790
##    y1[3]   15.11  15.14  2.59  2.19  10.78 19.36 1.01     2309     2993
##    y1[4]   19.06  19.02  2.52  2.09  14.90 23.28 1.01     3240     2807
##    y1[5]   22.94  22.91  2.54  2.07  18.85 27.28 1.00     3909     2335
##    y1[6]   26.83  26.79  2.61  2.11  22.52 31.14 1.01     3388     2375
##    y1[7]   30.71  30.67  2.54  2.12  26.43 34.93 1.01     2912     2913
##    y1[8]   34.73  34.68  2.62  2.24  30.39 39.07 1.00     2300     3063
##    y1[9]   38.59  38.57  2.65  2.36  34.25 43.03 1.00     1916     2091
##    y1[10]  42.44  42.46  2.79  2.52  38.01 47.09 1.01     1003     1584

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

outlier

ex10.stan

objective variable have outlier, y~cauchy(b0+b1*x,s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -16689 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17      -33.1588   0.000192288   0.000581333           1           1       47    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -33.16
##    b0         5.70
##    b1         1.93
##    s          3.18
##    m0[1]     11.48
##    m0[2]     19.72
##    m0[3]      9.90
##    m0[4]     10.86
##    m0[5]     22.74
##    m0[6]     15.07
##    m0[7]      7.14
##    m0[8]     19.06
##    m0[9]      7.24
##    m0[10]     7.93
##    m0[11]     6.06
##    m0[12]     6.90
##    m0[13]    12.03
##    m0[14]    18.47
##    m0[15]     9.32
##    m0[16]    11.59
##    m0[17]    16.29
##    m0[18]     7.16
##    m0[19]     6.45
##    m0[20]    20.68
##    y0[1]     15.93
##    y0[2]     18.77
##    y0[3]      6.73
##    y0[4]      9.68
##    y0[5]     24.91
##    y0[6]     16.19
##    y0[7]      8.19
##    y0[8]     19.66
##    y0[9]      3.80
##    y0[10]    11.99
##    y0[11]     3.38
##    y0[12]     6.69
##    y0[13]    10.58
##    y0[14]    22.80
##    y0[15]     7.16
##    y0[16]    10.62
##    y0[17]    16.49
##    y0[18]    11.78
##    y0[19]    10.87
##    y0[20]    18.60
##    m1[1]      6.06
##    m1[2]      7.92
##    m1[3]      9.77
##    m1[4]     11.62
##    m1[5]     13.48
##    m1[6]     15.33
##    m1[7]     17.18
##    m1[8]     19.04
##    m1[9]     20.89
##    m1[10]    22.74
##    y1[1]      8.61
##    y1[2]     11.55
##    y1[3]     12.65
##    y1[4]      9.79
##    y1[5]      7.93
##    y1[6]     15.34
##    y1[7]     19.99
##    y1[8]     18.29
##    y1[9]     18.94
##    y1[10]    27.20
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.60 -33.23 1.33 1.03 -36.20 -32.18 1.01      637      553
##    b0       5.71   5.61 1.35 1.28   3.67   8.06 1.02      263      282
##    b1       1.93   1.94 0.30 0.29   1.43   2.40 1.01      381      601
##    s        3.59   3.49 0.65 0.58   2.73   4.77 1.00     1478     1302
##    m0[1]   11.49  11.47 0.86 0.82  10.13  12.94 1.01      427      646
##    m0[2]   19.72  19.71 1.40 1.30  17.48  21.99 1.00     1611     1306
##    m0[3]    9.91   9.89 0.94 0.89   8.49  11.54 1.02      322      408
##    m0[4]   10.87  10.84 0.88 0.84   9.49  12.39 1.01      372      522
##    m0[5]   22.75  22.74 1.79 1.73  19.87  25.61 1.00     1058     1189
##    m0[6]   15.08  15.06 0.93 0.89  13.58  16.61 1.00     1739     1466
##    m0[7]    7.15   7.10 1.19 1.12   5.37   9.24 1.02      269      284
##    m0[8]   19.06  19.04 1.32 1.22  16.90  21.16 1.00     1724     1306
##    m0[9]    7.25   7.20 1.18 1.12   5.47   9.32 1.02      269      288
##    m0[10]   7.95   7.89 1.10 1.05   6.28   9.88 1.02      275      296
##    m0[11]   6.08   6.00 1.31 1.25   4.12   8.36 1.02      264      283
##    m0[12]   6.91   6.86 1.21 1.14   5.09   9.04 1.02      267      289
##    m0[13]  12.04  12.04 0.85 0.81  10.66  13.43 1.01      498      726
##    m0[14]  18.47  18.46 1.25 1.16  16.44  20.47 1.00     1828     1361
##    m0[15]   9.33   9.31 0.98 0.94   7.85  11.04 1.02      302      373
##    m0[16]  11.61  11.59 0.85 0.82  10.23  13.04 1.01      440      667
##    m0[17]  16.30  16.29 1.02 0.97  14.61  17.99 1.00     2007     1342
##    m0[18]   7.18   7.13 1.18 1.12   5.39   9.27 1.02      269      286
##    m0[19]   6.47   6.40 1.26 1.20   4.58   8.67 1.02      265      283
##    m0[20]  20.68  20.67 1.52 1.42  18.21  23.14 1.00     1432     1290
##    y0[1]   11.53  11.63 3.75 3.63   5.33  17.55 1.00     1518     1719
##    y0[2]   19.64  19.70 3.87 3.65  13.23  25.86 1.00     1929     1744
##    y0[3]   10.00  10.04 3.81 3.60   3.61  16.18 1.00     1629     1952
##    y0[4]   10.77  10.77 3.77 3.58   4.44  16.96 1.00     1615     1847
##    y0[5]   22.79  22.57 4.05 3.86  16.32  29.54 1.00     1976     2031
##    y0[6]   15.02  15.07 3.84 3.66   8.46  21.16 1.00     1931     1951
##    y0[7]    7.04   7.10 3.80 3.72   0.90  13.23 1.00     1387     1587
##    y0[8]   19.01  19.02 3.94 3.89  12.53  25.48 1.00     2043     1878
##    y0[9]    7.15   7.13 3.78 3.55   1.05  13.44 1.00     1012     1543
##    y0[10]   7.86   7.98 3.85 3.63   1.44  13.84 1.00     1512     1893
##    y0[11]   5.99   5.98 3.87 3.65  -0.30  12.26 1.00     1413     1824
##    y0[12]   6.91   6.95 3.96 3.86   0.49  13.39 1.00     1484     1534
##    y0[13]  11.96  11.95 3.69 3.68   5.96  17.82 1.00     1928     1779
##    y0[14]  18.30  18.28 3.91 3.69  11.96  24.92 1.00     2026     1719
##    y0[15]   9.26   9.30 3.83 3.83   2.94  15.56 1.00     1476     1770
##    y0[16]  11.49  11.58 3.76 3.53   5.35  17.34 1.00     2038     1899
##    y0[17]  16.32  16.39 3.79 3.65  10.16  22.53 1.00     1979     1922
##    y0[18]   7.31   7.26 3.81 3.59   1.07  13.59 1.00     1215     1536
##    y0[19]   6.50   6.52 3.89 3.77   0.22  12.71 1.00     1353     1758
##    y0[20]  20.59  20.44 3.92 3.84  14.22  27.11 1.00     1923     1772
##    m1[1]    6.08   6.00 1.31 1.25   4.12   8.36 1.02      264      283
##    m1[2]    7.93   7.88 1.10 1.05   6.26   9.87 1.02      275      296
##    m1[3]    9.78   9.76 0.94 0.90   8.35  11.43 1.02      317      393
##    m1[4]   11.64  11.62 0.85 0.82  10.27  13.07 1.01      444      667
##    m1[5]   13.49  13.48 0.85 0.80  12.11  14.88 1.00      882     1284
##    m1[6]   15.34  15.32 0.95 0.91  13.81  16.88 1.00     1821     1422
##    m1[7]   17.19  17.19 1.11 1.04  15.35  19.01 1.00     2007     1362
##    m1[8]   19.04  19.02 1.31 1.22  16.89  21.13 1.00     1728     1306
##    m1[9]   20.90  20.89 1.54 1.45  18.38  23.38 1.00     1373     1258
##    m1[10]  22.75  22.74 1.79 1.73  19.87  25.61 1.00     1058     1189
##    y1[1]    6.21   6.16 3.82 3.72   0.00  12.68 1.00     1157     1211
##    y1[2]    7.83   7.76 3.79 3.74   1.71  14.10 1.00     1471     1908
##    y1[3]    9.78   9.78 3.76 3.56   3.65  15.73 1.00     1489     1876
##    y1[4]   11.69  11.56 3.74 3.65   5.63  17.78 1.00     1873     1818
##    y1[5]   13.56  13.65 3.66 3.46   7.55  19.31 1.00     1792     1835
##    y1[6]   15.25  15.18 3.74 3.60   9.17  21.49 1.00     1931     2012
##    y1[7]   17.12  17.15 3.80 3.71  10.77  23.03 1.00     1814     1879
##    y1[8]   19.10  19.15 3.88 3.61  12.78  25.29 1.00     1807     1972
##    y1[9]   20.98  20.95 3.90 3.81  14.45  27.37 1.00     2027     1856
##    y1[10]  22.72  22.74 4.02 3.88  16.26  29.27 1.00     1900     1892

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -103.137 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       18      -9.24928   5.46446e-06   5.79417e-05     0.05057           1       33    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -9.25
##    b0         4.95
##    b1         1.99
##    s          0.51
##    m0[1]     10.91
##    m0[2]     19.41
##    m0[3]      9.29
##    m0[4]     10.27
##    m0[5]     22.52
##    m0[6]     14.62
##    m0[7]      6.44
##    m0[8]     18.73
##    m0[9]      6.54
##    m0[10]     7.26
##    m0[11]     5.33
##    m0[12]     6.19
##    m0[13]    11.48
##    m0[14]    18.12
##    m0[15]     8.68
##    m0[16]    11.03
##    m0[17]    15.87
##    m0[18]     6.46
##    m0[19]     5.73
##    m0[20]    20.39
##    y0[1]     12.35
##    y0[2]     14.07
##    y0[3]      9.38
##    y0[4]     10.42
##    y0[5]     22.60
##    y0[6]     15.03
##    y0[7]      6.64
##    y0[8]     18.72
##    y0[9]     14.66
##    y0[10]     6.47
##    y0[11]     5.70
##    y0[12]     6.02
##    y0[13]    11.99
##    y0[14]    18.18
##    y0[15]     8.40
##    y0[16]    11.08
##    y0[17]    15.59
##    y0[18]     5.81
##    y0[19]     8.51
##    y0[20]    16.12
##    m1[1]      5.33
##    m1[2]      7.24
##    m1[3]      9.15
##    m1[4]     11.06
##    m1[5]     12.97
##    m1[6]     14.88
##    m1[7]     16.79
##    m1[8]     18.70
##    m1[9]     20.61
##    m1[10]    22.52
##    y1[1]      0.65
##    y1[2]      7.29
##    y1[3]      9.15
##    y1[4]     12.71
##    y1[5]     12.55
##    y1[6]     15.05
##    y1[7]     17.45
##    y1[8]     18.93
##    y1[9]     20.80
##    y1[10]    22.20
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -11.50 -11.16   1.35 1.08 -14.17 -10.06 1.02      445      358
##    b0       4.95   4.96   0.31 0.31   4.44   5.42 1.01      897     1011
##    b1       1.98   1.98   0.07 0.06   1.86   2.08 1.01     1195     1397
##    s        0.61   0.60   0.19 0.17   0.36   0.96 1.01      735      284
##    m0[1]   10.88  10.88   0.21 0.22  10.52  11.21 1.01     1277     1309
##    m0[2]   19.32  19.35   0.34 0.30  18.72  19.85 1.01     2120     1554
##    m0[3]    9.26   9.26   0.22 0.23   8.89   9.61 1.01     1063     1208
##    m0[4]   10.24  10.25   0.21 0.22   9.88  10.58 1.01     1176     1277
##    m0[5]   22.42  22.45   0.43 0.39  21.65  23.10 1.01     1916     1669
##    m0[6]   14.56  14.57   0.24 0.22  14.15  14.93 1.00     2090     1526
##    m0[7]    6.43   6.44   0.27 0.28   5.97   6.86 1.01      916     1045
##    m0[8]   18.64  18.67   0.33 0.29  18.08  19.15 1.01     2162     1555
##    m0[9]    6.53   6.54   0.27 0.28   6.08   6.95 1.01      918     1069
##    m0[10]   7.24   7.25   0.26 0.26   6.82   7.64 1.01      938     1082
##    m0[11]   5.33   5.34   0.30 0.30   4.83   5.79 1.01      899     1012
##    m0[12]   6.18   6.19   0.28 0.28   5.72   6.61 1.01      910     1023
##    m0[13]  11.44  11.45   0.21 0.21  11.08  11.78 1.01     1386     1363
##    m0[14]  18.04  18.05   0.31 0.27  17.50  18.51 1.01     2197     1654
##    m0[15]   8.66   8.66   0.23 0.24   8.28   9.03 1.01     1009     1133
##    m0[16]  10.99  11.00   0.21 0.22  10.63  11.33 1.01     1297     1327
##    m0[17]  15.81  15.82   0.26 0.23  15.36  16.20 1.00     2219     1693
##    m0[18]   6.45   6.47   0.27 0.28   6.00   6.88 1.01      917     1045
##    m0[19]   5.73   5.74   0.29 0.30   5.25   6.18 1.01      905     1013
##    m0[20]  20.30  20.33   0.37 0.33  19.65  20.89 1.01     2050     1574
##    y0[1]   11.97  10.89  38.79 0.91   7.39  15.14 1.00     1960     2179
##    y0[2]   19.25  19.33  21.11 0.96  15.65  22.87 1.00     1938     1743
##    y0[3]    0.26   9.23 401.74 0.96   4.66  13.01 1.00     1887     1764
##    y0[4]    9.83  10.23  21.39 0.90   6.13  15.06 1.00     2039     1955
##    y0[5]   22.49  22.44  19.31 1.08  18.69  26.65 1.00     2171     1993
##    y0[6]   13.80  14.59 123.71 0.95  10.40  18.67 1.00     2069     2103
##    y0[7]    5.87   6.45  14.80 0.93   2.60  10.52 1.00     1925     1855
##    y0[8]   18.63  18.67   8.19 0.94  15.23  22.71 1.00     2011     2003
##    y0[9]    6.49   6.54  11.73 0.90   2.91  10.30 1.00     2171     1973
##    y0[10]   7.97   7.21  35.48 0.93   3.27  10.37 1.00     1975     1882
##    y0[11]   4.86   5.29  33.48 0.96   1.62   9.55 1.00     1836     1891
##    y0[12]   5.65   6.18  42.25 1.00   2.52  10.36 1.00     1932     1934
##    y0[13]  10.56  11.46  35.57 0.93   7.60  15.17 1.00     1981     1934
##    y0[14]  14.49  18.04 130.19 1.07  14.54  21.84 1.00     1974     1894
##    y0[15]  10.97   8.67  83.76 0.99   4.47  12.61 1.00     1949     1932
##    y0[16]  11.13  10.99  18.11 0.95   7.06  14.44 1.00     1680     1849
##    y0[17]   7.01  15.83 407.45 0.94  12.87  19.29 1.00     2018     2002
##    y0[18]   5.85   6.42  12.70 0.98   1.77   9.82 1.00     1848     1892
##    y0[19]   5.30   5.71  14.30 0.99   1.80   9.49 1.00     2019     1932
##    y0[20]  20.39  20.29  11.99 0.96  16.58  23.90 1.00     1958     1741
##    m1[1]    5.33   5.34   0.30 0.30   4.83   5.79 1.01      899     1012
##    m1[2]    7.23   7.24   0.26 0.27   6.80   7.63 1.01      937     1082
##    m1[3]    9.12   9.13   0.23 0.23   8.75   9.48 1.01     1050     1193
##    m1[4]   11.02  11.03   0.21 0.22  10.66  11.36 1.01     1303     1327
##    m1[5]   12.92  12.93   0.22 0.21  12.55  13.26 1.00     1742     1618
##    m1[6]   14.82  14.83   0.24 0.22  14.41  15.20 1.00     2129     1503
##    m1[7]   16.72  16.74   0.28 0.25  16.24  17.15 1.01     2243     1607
##    m1[8]   18.62  18.64   0.33 0.29  18.06  19.13 1.01     2164     1555
##    m1[9]   20.52  20.54   0.38 0.33  19.85  21.12 1.00     2035     1547
##    m1[10]  22.42  22.45   0.43 0.39  21.65  23.10 1.01     1916     1669
##    y1[1]    5.04   5.29  11.61 0.98   1.48   9.04 1.00     1774     1933
##    y1[2]    7.31   7.25  20.32 0.88   4.03  11.29 1.00     1663     1739
##    y1[3]    7.92   9.13  60.53 0.90   5.85  12.36 1.00     1831     1911
##    y1[4]    8.73  11.01  88.95 0.92   7.35  14.87 1.00     2065     2015
##    y1[5]   12.12  12.92  25.57 0.91   9.23  17.11 1.00     2228     2001
##    y1[6]   15.05  14.82  16.99 0.97  10.47  18.59 1.00     2053     1891
##    y1[7]   16.53  16.71  13.64 0.92  12.81  20.26 1.00     2026     1894
##    y1[8]   17.85  18.62  18.88 0.99  14.62  21.85 1.00     2071     1971
##    y1[9]   21.33  20.52  27.86 1.08  16.73  25.07 1.00     2033     1903
##    y1[10]  22.09  22.45  25.17 1.09  18.56  26.29 1.00     1993     1951

censored data

y i=1-N, y~N(m,s)  
  acutualy        ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11.stan

objective variable is censored

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan') 

mle=mdl$optimize(data=data)
## Initial log joint probability = -1957.23 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       30      -12.1858   0.000516082   1.11285e-06           1           1       42    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -12.19
##    b0        12.48
##    b1         2.52
##    s          2.78
##    m0[1]     20.66
##    m0[2]     21.02
##    m0[3]     18.89
##    m0[4]     23.78
##    m0[5]     22.97
##    m0[6]     30.78
##    m0[7]     22.05
##    m0[8]     23.67
##    y0[1]     20.27
##    y0[2]     22.04
##    y0[3]     20.26
##    y0[4]     24.80
##    y0[5]     26.36
##    y0[6]     27.19
##    y0[7]     19.73
##    y0[8]     23.66
##    m1[1]     12.60
##    m1[2]     15.06
##    m1[3]     17.52
##    m1[4]     19.98
##    m1[5]     22.45
##    m1[6]     24.91
##    m1[7]     27.37
##    m1[8]     29.83
##    m1[9]     32.30
##    m1[10]    34.76
##    y1[1]     15.40
##    y1[2]     15.01
##    y1[3]     14.48
##    y1[4]     18.98
##    y1[5]     20.79
##    y1[6]     26.78
##    y1[7]     22.70
##    y1[8]     29.94
##    y1[9]     29.30
##    y1[10]    35.01
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -13.14 -12.74 1.55 1.35 -16.09 -11.37 1.01      524      580
##    b0      12.34  12.45 5.11 4.73   3.67  20.88 1.03      206      230
##    b1       2.57   2.56 1.19 1.10   0.59   4.48 1.03      223      269
##    s        4.26   3.86 1.79 1.29   2.36   7.52 1.00      969     1052
##    m0[1]   20.66  20.65 1.89 1.68  17.57  23.70 1.02      327      499
##    m0[2]   21.03  21.00 1.81 1.60  18.08  23.91 1.01      376      615
##    m0[3]   18.86  18.89 2.43 2.17  14.76  22.84 1.02      234      299
##    m0[4]   23.84  23.80 1.69 1.48  21.19  26.53 1.00     1205     1053
##    m0[5]   23.01  22.98 1.62 1.38  20.45  25.53 1.00     1086     1007
##    m0[6]   30.95  31.01 4.14 3.74  24.42  37.56 1.02      324      427
##    m0[7]   22.07  22.05 1.65 1.39  19.38  24.66 1.01      649      879
##    m0[8]   23.72  23.69 1.68 1.46  21.07  26.34 1.00     1204     1077
##    y0[1]   20.59  20.60 4.86 4.09  12.44  28.11 1.01     1061     1568
##    y0[2]   21.16  21.06 4.63 4.09  13.99  28.65 1.00     1460     1825
##    y0[3]   18.78  18.81 5.23 4.39  10.03  26.86 1.00      809     1304
##    y0[4]   23.92  23.89 4.83 4.07  16.40  31.49 1.00     1948     1595
##    y0[5]   23.25  23.19 4.92 3.99  15.47  31.34 1.00     1922     1894
##    y0[6]   31.01  31.01 6.25 5.39  20.90  40.74 1.01      624     1092
##    y0[7]   22.11  22.17 4.92 4.22  14.61  29.60 1.00     1743     1693
##    y0[8]   23.68  23.66 4.92 4.01  15.58  31.23 1.00     2221     1543
##    m1[1]   12.46  12.56 5.06 4.67   3.86  20.92 1.03      206      224
##    m1[2]   14.97  15.06 3.97 3.61   8.03  21.45 1.03      207      227
##    m1[3]   17.47  17.51 2.94 2.64  12.47  22.33 1.03      215      274
##    m1[4]   19.97  19.97 2.07 1.78  16.60  23.33 1.02      273      399
##    m1[5]   22.48  22.43 1.62 1.36  19.89  24.99 1.00      874      919
##    m1[6]   24.98  24.94 1.92 1.66  21.88  27.96 1.00     1015      916
##    m1[7]   27.49  27.49 2.73 2.42  23.09  31.70 1.01      484      669
##    m1[8]   29.99  30.03 3.73 3.38  24.16  35.91 1.02      345      441
##    m1[9]   32.49  32.53 4.81 4.37  24.91  40.02 1.02      300      402
##    m1[10]  35.00  35.05 5.92 5.45  25.52  44.30 1.02      278      377
##    y1[1]   12.41  12.39 6.93 6.03   1.73  23.03 1.02      343      684
##    y1[2]   14.92  14.92 6.22 5.21   4.98  24.95 1.01      389      866
##    y1[3]   17.62  17.52 5.72 4.70   9.03  26.87 1.01      510      792
##    y1[4]   20.01  19.89 5.25 4.25  11.62  28.00 1.00     1043     1541
##    y1[5]   22.69  22.56 4.96 4.23  15.07  30.68 1.00     1586     1770
##    y1[6]   24.72  24.82 4.87 4.01  17.00  32.45 1.00     1877     1563
##    y1[7]   27.31  27.25 5.37 4.49  19.13  35.86 1.00     1276     1398
##    y1[8]   30.00  29.91 5.95 5.05  20.62  39.30 1.01      660     1403
##    y1[9]   32.48  32.42 6.78 5.65  21.31  43.36 1.01      474      846
##    y1[10]  35.04  35.19 7.50 6.36  22.83  47.20 1.01      429      857

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -177.751 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       27      -19.1875   0.000483383   0.000105904      0.9521      0.9521       33    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -19.19
##    b0         5.85
##    b1         4.16
##    s          4.27
##    m0[1]     19.32
##    m0[2]     19.92
##    m0[3]     16.40
##    m0[4]     24.47
##    m0[5]     23.13
##    m0[6]     35.99
##    m0[7]     21.61
##    m0[8]     24.28
##    y0[1]     16.55
##    y0[2]     28.06
##    y0[3]      9.68
##    y0[4]     27.35
##    y0[5]     26.15
##    y0[6]     32.04
##    y0[7]     18.23
##    y0[8]     22.68
##    m1[1]      6.04
##    m1[2]     10.09
##    m1[3]     14.15
##    m1[4]     18.21
##    m1[5]     22.27
##    m1[6]     26.32
##    m1[7]     30.38
##    m1[8]     34.44
##    m1[9]     38.49
##    m1[10]    42.55
##    y1[1]      7.09
##    y1[2]      8.03
##    y1[3]     14.02
##    y1[4]     20.92
##    y1[5]     22.69
##    y1[6]     29.69
##    y1[7]     28.38
##    y1[8]     33.89
##    y1[9]     36.67
##    y1[10]    41.44
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median    sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -19.63 -19.23  1.55 1.24 -22.78 -17.98 1.01      295      264
##    b0       1.58   2.59  7.05 6.08 -11.64  10.52 1.04      182      146
##    b1       5.23   4.92  1.63 1.32   3.24   8.29 1.03      197      149
##    s        6.55   5.82  2.92 2.06   3.57  11.86 1.01      368      345
##    m0[1]   18.54  18.77  2.52 2.18  14.23  22.13 1.03      282      225
##    m0[2]   19.28  19.47  2.40 2.09  15.26  22.82 1.03      314      250
##    m0[3]   14.86  15.25  3.32 2.80   8.86  19.42 1.04      202      175
##    m0[4]   25.01  24.81  2.13 1.72  22.09  28.57 1.01     1076      796
##    m0[5]   23.33  23.27  2.06 1.74  20.29  26.83 1.01      749      743
##    m0[6]   39.51  38.57  5.47 4.28  32.96  49.49 1.01      274      270
##    m0[7]   21.41  21.45  2.13 1.82  18.06  24.87 1.02      459      430
##    m0[8]   24.77  24.58  2.11 1.70  21.84  28.33 1.01     1040      751
##    y0[1]   18.55  18.54  7.75 6.20   6.89  29.77 1.00     1778     1370
##    y0[2]   19.36  19.63  7.43 6.41   7.37  30.57 1.00     1275     1216
##    y0[3]   14.77  15.24  7.68 6.58   2.26  25.87 1.01      965      872
##    y0[4]   24.92  24.84  7.66 6.03  12.68  36.96 1.00     1892     1523
##    y0[5]   23.26  23.29  7.27 6.06  11.81  34.60 1.00     1783     1205
##    y0[6]   39.66  38.58  9.22 7.34  27.50  55.61 1.01      624      469
##    y0[7]   21.08  21.36  7.70 6.10   9.55  32.49 1.00     1564     1097
##    y0[8]   24.89  24.66  7.54 6.15  13.55  36.76 1.00     2049     1618
##    m1[1]    1.82   2.83  6.98 6.02 -11.28  10.68 1.04      182      146
##    m1[2]    6.93   7.71  5.48 4.60  -3.22  13.99 1.04      184      152
##    m1[3]   12.03  12.53  4.05 3.46   4.70  17.37 1.04      189      164
##    m1[4]   17.13  17.42  2.80 2.37  12.26  21.05 1.04      232      208
##    m1[5]   22.24  22.22  2.08 1.77  19.06  25.71 1.02      556      612
##    m1[6]   27.34  27.05  2.42 1.93  24.21  31.55 1.01      910      695
##    m1[7]   32.44  31.91  3.53 2.79  28.06  38.91 1.01      397      323
##    m1[8]   37.55  36.76  4.91 3.85  31.62  46.59 1.01      292      272
##    m1[9]   42.65  41.67  6.39 5.05  34.92  54.45 1.02      256      276
##    m1[10]  47.75  46.50  7.92 6.35  38.15  62.62 1.02      239      233
##    y1[1]    2.01   3.23  9.52 8.07 -14.36  15.15 1.02      386      349
##    y1[2]    6.71   7.86  9.19 7.78  -9.60  18.74 1.01      469      405
##    y1[3]   12.02  12.48  8.10 6.74  -2.01  23.75 1.01      624     1014
##    y1[4]   17.05  17.48  8.00 6.21   3.27  29.08 1.01      871      775
##    y1[5]   22.13  22.20  7.48 6.11   9.89  34.17 1.00     1842     1544
##    y1[6]   27.43  26.97  7.51 6.23  16.14  40.22 1.00     1730     1061
##    y1[7]   32.48  32.04  7.85 6.32  20.76  45.86 1.00     1370      722
##    y1[8]   37.75  36.76  8.87 6.96  25.99  52.87 1.00      936      708
##    y1[9]   43.00  41.74  9.45 7.58  30.54  60.08 1.00      527      527
##    y1[10]  47.63  46.20 10.60 8.59  34.35  66.00 1.01      419      507

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

sensitivity/specificity analysis

ex14.stan

estimating sens and spec

data {
  int N;
  array[N] int x;
  array[N] int y;
}
parameters {
  real<lower=0,upper=1> p;
  real<lower=0,upper=1> se;
  real<lower=0,upper=1> sp;
}
model {
  p~uniform(0,1);
  se~uniform(0,1);
  sp~uniform(0,1);
  for (i in 1:N) {
    y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
  }
}
generated quantities {
  real ppv;
  real npv;
  ppv=se*p/((se*p)+((1-p)*(1-sp)));
  npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}
n=20
data=list(N=n,
          x=sample(0:1,n,replace=T),
          y=sample(0:1,n,replace=T))

mdl=cmdstan_model('./ex14.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -24.7567 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5      -13.4602    0.00123284   6.67477e-06           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -13.46
##      p        0.32
##      se       0.60
##      sp       0.40
##      ppv      0.32
##      npv      0.68
system.time({
  mcmc=goMCMC(mdl,data)
  seeMCMC(mcmc,ch=F)
})
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -19.41 -19.11 1.34 1.14 -22.09 -17.89 1.00      672      917
##      p      0.51   0.50 0.30 0.39   0.04   0.96 1.00      543      386
##      se     0.58   0.58 0.14 0.15   0.34   0.80 1.00     2833     1602
##      sp     0.42   0.41 0.14 0.15   0.20   0.65 1.00     2813     1509
##      ppv    0.51   0.51 0.30 0.39   0.05   0.96 1.01      567      410
##      npv    0.49   0.49 0.30 0.41   0.04   0.95 1.00      579      448

##   ユーザ システム     経過 
##    0.557    0.116    0.695

appendix: bimodal distribution, mixed normal distribution

stan

data {
  int<lower=0> N;
  real y[N];
}
parameters {
  real<lower=0, upper=1> theta; //ratio of mix
  real mu1;
  real mu2;
  real<lower=0> sigma1;
  real<lower=0> sigma2;
}
model {
  for (n in 1:N) {
    target += log_mix(theta,
                      normal_lpdf(y[n] | mu1, sigma1),
                      normal_lpdf(y[n] | mu2, sigma2));
  }
}

EM algorithm

y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 3 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -250 100  6 -528 -530
## 
## Clustering table:
##  1  2  3 
## 33 25 42
rst$parameters
## $pro
## [1] 0.329 0.254 0.417
## 
## $mean
##      1      2      3 
## -5.336  0.201  5.245 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 3
## 
## $variance$sigmasq
## [1] 1.05
plot(rst)